##############################################################

Tumor=$1
Normal=$2
line=$3
config_file=$4

source ${config_file}

##############################################################
## 定义bam文件以及输出的文件夹

tumor_bam=${vcf_path}/${Tumor}_${Normal}.mutect2.bamout.bam
vcf_qc_tmp_path=${vcf_qc_diff_path}/tmp/${Tumor}_${Normal}
vcf_qc_path=${vcf_qc_diff_path}

##############################################################
## line的格式：chr1:1202552:C:T

pos=`echo ${line} | awk -F':' '{print $1":"$2"-"$2}' `
ref=`echo ${line} | awk -F':' '{print $3}' `
alt=`echo ${line} | awk -F':' '{print $4}' `
start_pos=`echo ${line} | awk -F':' '{print $2}' `

ref_length=`expr length $ref`
alt_length=`expr length $alt`

##############################################################
#### 1、判断突变情况，插入、缺失、单碱基的替换

if [ ${ref_length} -gt ${alt_length} ]
then
## 缺失
mut_type="D"
mut_length=` expr ${ref_length} - ${alt_length} `
mut_type_forDistance=${mut_length}${mut_type}
elif [ ${ref_length} -lt ${alt_length} ]
then
## 插入
mut_type="I"
mut_length=` expr ${alt_length} - ${ref_length} `
mut_type_forDistance=${mut_length}${mut_type}
elif [ ${ref_length} -eq ${alt_length} ]
then
## SNV的改变
if [ ${ref_length} -eq 1 ]
then
mut_type="SNV"
mut_length=` expr ${alt_length} - ${ref_length} `
mut_base_pos=${start_pos}
mut_type_forDistance=${mut_type}
## 多碱基的替换
elif [ ${ref_length} -gt 1 ]
then
mut_type="MBS"
mut_length=${ref_length}
mut_type_forDistance=${mut_type}
fi

fi


##############################################################
#### 2、提取该位置的所有reads
## https://gatk.broadinstitute.org/hc/en-us/articles/360037593851-Mutect2
## 依据mutect2去call突变前的质控，去除不合格的reads
## 1、去除mutect2组装
## 2、去除比对到多条染色体的reads（Mutect2未去除）
## 3、只考虑两个配对的reads都匹配上的
## 	F 4:该read没有比对上 ; 8:与该read配对的mate read没有比对上
## 4、去除重复reads
## 	-F 1024
## 	https://www.biostars.org/p/208897/
## 	https://gatk.broadinstitute.org/hc/en-us/articles/360037592051-NotDuplicateReadFilter
## 5、去除NotSecondaryAlignmentReadFilter
## 	-F 256
## 	https://gatk.broadinstitute.org/hc/en-us/articles/360037227412-NotSecondaryAlignmentReadFilter
## 6、去除PassesVendorQualityCheckReadFilter
## 	-F 512
##  Filter out reads failing platform/vendor quality checks 
## 	https://gatk.broadinstitute.org/hc/en-us/articles/360037226972-PassesVendorQualityCheckReadFilter
## 7、去除
##	-F 2048
## 	https://www.biostars.org/p/8318/

${samtools} view -F 12 -F 1024 -F 256 -F 512 -F 2048 -h ${tumor_bam} ${pos} | grep -v ArtificialHaplotypeRG | \
awk -F'\t' '{OFS="\t"}{
if($0~"^@"){print};
if($0!~"^@" && $7=="="){print}}' | \
${samtools} view -b -h - > ${vcf_qc_tmp_path}/${Tumor}_${Normal}.${pos}.bam
${samtools} index ${vcf_qc_tmp_path}/${Tumor}_${Normal}.${pos}.bam


##############################################################
#### 3、得到突变reads的name
## 1、对于INDEL，因为Sam2Tsv会按照每个位点拆分bam文件，通过同一reads的连续位点的插入缺失的突变类型以及突变的长度，确认我们感兴趣的突变
## 2、对于SNV，按照突变位置以及碱基的替换之间去匹配
## 3、对MBS，仿照INDEL的思路
## 4、reads同时包含了flag的信息，以保证完全匹配(存在相同reads，正义链和反义链只有一条有突变，保留有突变的reads)

## http://lindenb.github.io/jvarkit/Sam2Tsv.html

###################################
## INS的确认
if [ ${mut_type} == "I" ]
then
java -jar ${tools_path}/jvarkit/dist/sam2tsv.jar -R ${ref_fasta} ${vcf_qc_tmp_path}/${Tumor}_${Normal}.${pos}.bam \
> ${vcf_qc_tmp_path}/${Tumor}_${Normal}.sam2tsv.${pos}.all.list

## 1、提取突变位置开始，所有reads后面丢失长度相同的碱基
# cat ${vcf_qc_tmp_path}/${Tumor}_${Normal}.sam2tsv.${pos}.all.list | grep -A ${mut_length} ${start_pos} | awk -F'\t' '{print $6"\t"$1":"$2}'\
# > ${vcf_qc_tmp_path}/${Tumor}_${Normal}.sam2tsv.${pos}.altpos.list

## 2、当突变的长度过长，非所有reads都能满足该长度，导致开始位置往后该长度到其它reads，会引起错误，使用下面的版本
## 	A、每个reads，计算该位置到突变的位置长度，不满足突变长度则到突变的最后位置
## 	B、提取包含插入的reads，提高速度
## 	C、reads需要满足，突变的reads一致；突变类型一致；以访该区域存在原先插入的reads；如原先为该区域为CAA，然后插入一个A变成CAAA，突变碱基捕获则认为和原来一样，已解决；

echo "--:" > ${vcf_qc_tmp_path}/${Tumor}_${Normal}.sam2tsv.${pos}.altpos.list

for read in `cat ${vcf_qc_tmp_path}/${Tumor}_${Normal}.sam2tsv.${pos}.all.list | awk -F'\t' '{if($10==mut_type)print $1":"$2}' mut_type=${mut_type} | uniq`
do

cat ${vcf_qc_tmp_path}/${Tumor}_${Normal}.sam2tsv.${pos}.all.list | \
awk -F'\t' '{match_id=$1":"$2; if(match_id ~ read )print}'  read=${read} | \
grep -A ${mut_length} ${start_pos} | \
awk -F'\t' '{print $6"\t"$1":"$2"\t"$10}' \
>> ${vcf_qc_tmp_path}/${Tumor}_${Normal}.sam2tsv.${pos}.altpos.list
echo "--:" >> ${vcf_qc_tmp_path}/${Tumor}_${Normal}.sam2tsv.${pos}.altpos.list

done

## 每个reads，匹配相同起始位置且突变完全相同的碱基
## 产生插入的reads的标签
indel_base=$(printf "%-${mut_length}s" "I")
indel_base=`echo "${indel_base// /I}"`
indel_base="M"${indel_base}

alt_readname=""
for read in `cat ${vcf_qc_tmp_path}/${Tumor}_${Normal}.sam2tsv.${pos}.altpos.list | awk -F'\t' '{print $2}' | grep -v "\-\-\:" | uniq`
do
tmp=`cat ${vcf_qc_tmp_path}/${Tumor}_${Normal}.sam2tsv.${pos}.altpos.list | grep ${read} | awk -F'\t' '{print $1}'`
tmp=`echo $tmp | sed 's/ //g'`

tmp_2=`cat ${vcf_qc_tmp_path}/${Tumor}_${Normal}.sam2tsv.${pos}.altpos.list | grep ${read} | awk -F'\t' '{print $3}'`
tmp_2=`echo $tmp_2 | sed 's/ //g'`

if [ ${tmp} == ${alt} ] && [ ${tmp_2} == ${indel_base} ]
then
alt_readname=${read}"|"${alt_readname}
fi
done
alt_readname=`echo ${alt_readname} | sed 's/|$//'`
fi

###################################
## DEL的确认
if [ ${mut_type} == "D" ]
then
java -jar ${tools_path}/jvarkit/dist/sam2tsv.jar -R ${ref_fasta} ${vcf_qc_tmp_path}/${Tumor}_${Normal}.${pos}.bam \
> ${vcf_qc_tmp_path}/${Tumor}_${Normal}.sam2tsv.${pos}.all.list

## 1、提取突变位置开始，所有reads后面丢失长度相同的碱基
# cat ${vcf_qc_tmp_path}/${Tumor}_${Normal}.sam2tsv.${pos}.all.list | grep -A ${mut_length} ${start_pos} | awk -F'\t' '{print $6"\t"$1":"$2}'\
# > ${vcf_qc_tmp_path}/${Tumor}_${Normal}.sam2tsv.${pos}.altpos.list

## 2、当突变的长度过长，非所有reads都能满足该长度，导致开始位置往后该长度到其它reads，会引起错误，使用第二版本
## 每个reads，计算该位置到最后位置的长度
echo "--:" > ${vcf_qc_tmp_path}/${Tumor}_${Normal}.sam2tsv.${pos}.altpos.list

for read in `cat ${vcf_qc_tmp_path}/${Tumor}_${Normal}.sam2tsv.${pos}.all.list | awk -F'\t' '{if($10==mut_type)print $1":"$2}' mut_type=${mut_type} | uniq`
do

cat ${vcf_qc_tmp_path}/${Tumor}_${Normal}.sam2tsv.${pos}.all.list | \
awk -F'\t' '{match_id=$1":"$2; if(match_id ~ read )print}'  read=${read} | \
grep -A ${mut_length} ${start_pos} | \
awk -F'\t' '{print $6"\t"$1":"$2}' \
>> ${vcf_qc_tmp_path}/${Tumor}_${Normal}.sam2tsv.${pos}.altpos.list
echo "--:" >> ${vcf_qc_tmp_path}/${Tumor}_${Normal}.sam2tsv.${pos}.altpos.list

done

## 产生缺失的reads
del_base=$(printf "%-${mut_length}s" ".")
del_base=`echo "${del_base// /.}"`
del_base=${alt}${del_base}

alt_readname=""
for read in `cat ${vcf_qc_tmp_path}/${Tumor}_${Normal}.sam2tsv.${pos}.altpos.list | awk -F'\t' '{print $2}' | grep -v "\-\-\:" | uniq`
do
tmp=`cat ${vcf_qc_tmp_path}/${Tumor}_${Normal}.sam2tsv.${pos}.altpos.list | grep ${read} | awk -F'\t' '{print $1}'`
tmp=`echo $tmp | sed 's/ //g'`
if [ ${tmp} == ${del_base} ]
then
alt_readname=${read}"|"${alt_readname}
fi
done

alt_readname=`echo ${alt_readname} | sed 's/|$//'`
fi

###################################
## SNV的确认
if [ ${mut_type} == "SNV" ]
then
java -jar ${tools_path}/jvarkit/dist/sam2tsv.jar -R ${ref_fasta} ${vcf_qc_tmp_path}/${Tumor}_${Normal}.${pos}.bam | \
awk -F'\t' '{if($8==mut_base_pos && $6==alt )print}' mut_type=${mut_type} mut_base_pos=${mut_base_pos} alt=${alt} \
> ${vcf_qc_tmp_path}/${Tumor}_${Normal}.sam2tsv.${pos}.list
alt_readname=`cat ${vcf_qc_tmp_path}/${Tumor}_${Normal}.sam2tsv.${pos}.list | awk -F '\t' '{print $1":"$2}' | uniq | tr '\n' '|' | sed 's/|$//' `
fi

###################################
## MBS
## DBS以及三个及以上的变化
if [ ${mut_type} == "MBS" ]
then

java -jar ${tools_path}/jvarkit/dist/sam2tsv.jar -R ${ref_fasta} ${vcf_qc_tmp_path}/${Tumor}_${Normal}.${pos}.bam \
> ${vcf_qc_tmp_path}/${Tumor}_${Normal}.sam2tsv.${pos}.all.list

## 提取突变位置开始，所有reads后面插入大小的位置上的碱基

match_length=`expr ${mut_length} - 1`

echo "--:" > ${vcf_qc_tmp_path}/${Tumor}_${Normal}.sam2tsv.${pos}.altpos.list
for read in `cat ${vcf_qc_tmp_path}/${Tumor}_${Normal}.sam2tsv.${pos}.all.list | awk -F'\t' '{if($10==mut_type)print $1":"$2}' mut_type="M" | uniq`
do
cat ${vcf_qc_tmp_path}/${Tumor}_${Normal}.sam2tsv.${pos}.all.list | \
awk -F'\t' '{match_id=$1":"$2; if(match_id ~ read )print}'  read=${read} | \
grep -A ${match_length} ${start_pos} | \
awk -F'\t' '{print $6"\t"$1":"$2}' \
>> ${vcf_qc_tmp_path}/${Tumor}_${Normal}.sam2tsv.${pos}.altpos.list
echo "--:" >> ${vcf_qc_tmp_path}/${Tumor}_${Normal}.sam2tsv.${pos}.altpos.list
done

## 每个reads，匹配相同起始位置且突变完全相同的碱基
alt_readname=""
for read in `cat ${vcf_qc_tmp_path}/${Tumor}_${Normal}.sam2tsv.${pos}.altpos.list | awk -F'\t' '{print $2}' | grep -v "\-\-\:" | uniq`
do
tmp=`cat ${vcf_qc_tmp_path}/${Tumor}_${Normal}.sam2tsv.${pos}.altpos.list | grep ${read} | awk -F'\t' '{print $1}'`
tmp=`echo $tmp | sed 's/ //g'`
if [ ${tmp} == ${alt} ]
then
alt_readname=${read}"|"${alt_readname}
fi
done
alt_readname=`echo ${alt_readname} | sed 's/|$//'`

fi


##############################################################
## 若突变所在的所有reads，不满足质控的条件
if [ -z ${alt_readname} ]
then
echo "No Alt Reads"
echo "${line},NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA" >> ${vcf_qc_path}/${Tumor}_${Normal}_MutQc.csv
elif [ -n ${alt_readname} ]
then

## 判读突变reads在多少tumor中
Tumor_alt_counts=`${samtools} view ${vcf_qc_tmp_path}/${Tumor}_${Normal}.${pos}.bam | \
awk -F'\t' '{match_id=$1":"$2; if(match_id ~ alt_readname )print}' alt_readname=${alt_readname} | \
grep "RG:Z:${Tumor}" | wc -l`

## 若Tumor中无满足质控标准的突变的reads
if [ ${Tumor_alt_counts} -eq 0 ]
then
echo "${line},NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA" >> ${vcf_qc_path}/${Tumor}_${Normal}_MutQc.csv
else

##############################################################
####  4、reads的数量
# 总数量
tumor_counts=`${samtools} view ${vcf_qc_tmp_path}/${Tumor}_${Normal}.${pos}.bam | grep "RG:Z:${Tumor}" | wc -l`
normal_counts=`${samtools} view ${vcf_qc_tmp_path}/${Tumor}_${Normal}.${pos}.bam | grep "RG:Z:${Normal}" | wc -l`

## 正义链和反义链的reads
## https://broadinstitute.github.io/picard/explain-flags.html
## https://www.biostars.org/p/14378/
## ${samtools} view -F 20 ... : forward strand( read unmapped ,  read reverse strand)
# ${samtools} view -f 16 ... : reverse strand( read reverse strand)
# ${samtools} view -f 4 ... : unmapped
Tumor_read1_counts=`${samtools} view -f 16 ${vcf_qc_tmp_path}/${Tumor}_${Normal}.${pos}.bam | grep "RG:Z:${Tumor}" | wc -l`
Tumor_read2_counts=`${samtools} view -F 20  ${vcf_qc_tmp_path}/${Tumor}_${Normal}.${pos}.bam | grep "RG:Z:${Tumor}" | wc -l`


##############################################################
####  5、提取tumor的alt的reads、normal的alt的reads、tumor的ref的reads、normal的ref的reads

##############################
## 1、提取Tumor突变的reasds
## 头行
${samtools} view -H ${vcf_qc_tmp_path}/${Tumor}_${Normal}.${pos}.bam > ${vcf_qc_tmp_path}/${Tumor}_${Normal}.${pos}.altreads.tumor.sam

## 提取Tumor中突变的reads
${samtools} view ${vcf_qc_tmp_path}/${Tumor}_${Normal}.${pos}.bam | \
awk -F'\t' '{match_id=$1":"$2; if(match_id ~ alt_readname )print}'  alt_readname=${alt_readname} | \
grep "RG:Z:${Tumor}" \
>> ${vcf_qc_tmp_path}/${Tumor}_${Normal}.${pos}.altreads.tumor.sam

${samtools} view -b -h ${vcf_qc_tmp_path}/${Tumor}_${Normal}.${pos}.altreads.tumor.sam > ${vcf_qc_tmp_path}/${Tumor}_${Normal}.${pos}.altreads.tumor.bam
${samtools} index ${vcf_qc_tmp_path}/${Tumor}_${Normal}.${pos}.altreads.tumor.bam

##############################
## 2、提取Numor突变的reasds
## 头行
${samtools} view -H ${vcf_qc_tmp_path}/${Tumor}_${Normal}.${pos}.bam > ${vcf_qc_tmp_path}/${Tumor}_${Normal}.${pos}.altreads.normal.sam

## 提取Normal中突变的reads
${samtools} view ${vcf_qc_tmp_path}/${Tumor}_${Normal}.${pos}.bam | \
awk -F'\t' '{match_id=$1":"$2; if(match_id ~ alt_readname )print}'  alt_readname=${alt_readname} | \
grep "RG:Z:${Normal}" \
>> ${vcf_qc_tmp_path}/${Tumor}_${Normal}.${pos}.altreads.normal.sam

${samtools} view -b -h ${vcf_qc_tmp_path}/${Tumor}_${Normal}.${pos}.altreads.normal.sam > ${vcf_qc_tmp_path}/${Tumor}_${Normal}.${pos}.altreads.normal.bam
${samtools} index ${vcf_qc_tmp_path}/${Tumor}_${Normal}.${pos}.altreads.normal.bam

##############################
## 3、提取Tumor的参考的reasds
## 头行
${samtools} view -H ${vcf_qc_tmp_path}/${Tumor}_${Normal}.${pos}.bam > ${vcf_qc_tmp_path}/${Tumor}_${Normal}.${pos}.refreads.tumor.sam

## 去除突变的reads
${samtools} view ${vcf_qc_tmp_path}/${Tumor}_${Normal}.${pos}.bam | \
awk -F'\t' '{match_id=$1":"$2; if(match_id !~ alt_readname )print}'  alt_readname=${alt_readname} | \
grep "RG:Z:${Tumor}" \
>> ${vcf_qc_tmp_path}/${Tumor}_${Normal}.${pos}.refreads.tumor.sam

## 转化为bam文件
${samtools} view -b -h ${vcf_qc_tmp_path}/${Tumor}_${Normal}.${pos}.refreads.tumor.sam > ${vcf_qc_tmp_path}/${Tumor}_${Normal}.${pos}.refreads.tumor.bam
${samtools} index ${vcf_qc_tmp_path}/${Tumor}_${Normal}.${pos}.refreads.tumor.bam

##############################
## 4、提取Numor的参考的reasds
## 头行
${samtools} view -H ${vcf_qc_tmp_path}/${Tumor}_${Normal}.${pos}.bam > ${vcf_qc_tmp_path}/${Tumor}_${Normal}.${pos}.refreads.normal.sam

## 去除突变的reads
${samtools} view ${vcf_qc_tmp_path}/${Tumor}_${Normal}.${pos}.bam | \
awk -F'\t' '{match_id=$1":"$2; if(match_id !~ alt_readname )print}'  alt_readname=${alt_readname} | \
grep "RG:Z:${Normal}" \
>> ${vcf_qc_tmp_path}/${Tumor}_${Normal}.${pos}.refreads.normal.sam

## 转化为bam文件
${samtools} view -b -h ${vcf_qc_tmp_path}/${Tumor}_${Normal}.${pos}.refreads.normal.sam > ${vcf_qc_tmp_path}/${Tumor}_${Normal}.${pos}.refreads.normal.bam
${samtools} index ${vcf_qc_tmp_path}/${Tumor}_${Normal}.${pos}.refreads.normal.bam

##############################################################
## 6、判断突变在各个链的reads数量
Tumor_read1_alt_counts=`${samtools} view -f 16 ${vcf_qc_tmp_path}/${Tumor}_${Normal}.${pos}.altreads.tumor.bam | wc -l`
Tumor_read2_alt_counts=`${samtools} view -F 20 ${vcf_qc_tmp_path}/${Tumor}_${Normal}.${pos}.altreads.tumor.bam | wc -l`


##############################################################
## 7、计算每个突变距离最近端的距离，每个突变距离最近端的中位距离的绝对偏差的中位数
# DistanceToAlignmentEndMedian
# DistanceToAlignmentEndMAD
${samtools} view ${vcf_qc_tmp_path}/${Tumor}_${Normal}.${pos}.altreads.tumor.bam | awk -F'\t' '{print $6}' \
> ${vcf_qc_tmp_path}/${Tumor}_${Normal}.${pos}.distance.txt

## 提取包含感兴趣的突变reads的mpileup
## -A, --count-orphans
## Do not skip anomalous read pairs in variant calling. 
## Anomalous read pairs are those marked in the FLAG field as paired in sequencing but without the properly-paired flag set.

## BAQ (Base Alignment Quality)
## BAQ is the Phred-scaled probability of a read base being misaligned. It greatly helps to reduce false SNPs caused by misalignments. 
## BAQ is calculated using the probabilistic realignment method described in the paper “Improving SNP discovery by base alignment quality”, Heng Li, Bioinformatics, Volume 27, Issue 8 
## <https://doi.org/10.1093/bioinformatics/btr076>
## BAQ is turned on when a reference file is supplied using the -f option. To disable it, use the -B option.
${samtools} mpileup --min-BQ 0 --min-MQ 0 -B -d10000 --count-orphans -sf ${ref_fasta} --output-BP -r ${pos} ${vcf_qc_tmp_path}/${Tumor}_${Normal}.${pos}.altreads.tumor.bam \
> ${vcf_qc_tmp_path}/${Tumor}_${Normal}.mpileup.${pos}.list 


## 产生reads的突变位置
cat ${vcf_qc_tmp_path}/${Tumor}_${Normal}.mpileup.${pos}.list | awk -F'\t' '{print $8}' | tr ',' '\n' \
> ${vcf_qc_tmp_path}/${Tumor}_${Normal}.${pos}.mutpos.txt

## 插入和缺失，CIGAR包含信息，按照插入缺失的位置计算
## 单碱基改变，CIGAR无信息，按照突变出现在reads的位置
## MBS，CIGAR无信息，按照突变出现在reads的位置，右端减去突变的长度
${Rscript} ${scripts_path}/DistanceToAlignmentEnd.R \
--input_file ${vcf_qc_tmp_path}/${Tumor}_${Normal}.${pos}.distance.txt \
--mutpos_file ${vcf_qc_tmp_path}/${Tumor}_${Normal}.${pos}.mutpos.txt \
--out_file ${vcf_qc_tmp_path}/${Tumor}_${Normal}.${pos}.distance.calculate.txt \
--mutType ${mut_type_forDistance} \
--mut_length ${mut_length}

DistanceToAlignmentEndMedian=`cat ${vcf_qc_tmp_path}/${Tumor}_${Normal}.${pos}.distance.calculate.txt | sed '1d' | awk -F'\t' '{print $1}'`
DistanceToAlignmentEndMAD=`cat ${vcf_qc_tmp_path}/${Tumor}_${Normal}.${pos}.distance.calculate.txt | sed '1d' | awk -F'\t' '{print $2}'`

##############################################################
# 8、minmapqualitydifference
## Maximum difference between median mapping quality of ALT reads in tumor and REF reads in normal
medianmapquality_altTumor=`${samtools} view ${vcf_qc_tmp_path}/${Tumor}_${Normal}.${pos}.altreads.tumor.bam | awk -F'\t' '{print $5}' | ${datamash} median 1`
medianmapquality_refNormal=`${samtools} view ${vcf_qc_tmp_path}/${Tumor}_${Normal}.${pos}.refreads.normal.bam | awk -F'\t' '{print $5}' | ${datamash} median 1`

## expr 无法计算非整数，使用awk，取绝对值
minmapqualitydifference=`echo "" | \
awk '{ 
minmapqualitydifference = (medianmapquality_altTumor - medianmapquality_refNormal);
if(minmapqualitydifference >= 0){minmapqualitydifference=minmapqualitydifference};
if(minmapqualitydifference < 0){minmapqualitydifference=-minmapqualitydifference};
print minmapqualitydifference
}' \
medianmapquality_altTumor=${medianmapquality_altTumor} medianmapquality_refNormal=${medianmapquality_refNormal}`

medianmapquality_Tumor=`${samtools} view ${vcf_qc_tmp_path}/${Tumor}_${Normal}.${pos}.bam | grep "RG:Z:${Tumor}" | awk -F'\t' '{print $5}' | ${datamash} median 1`
medianmapquality_Normal=`${samtools} view ${vcf_qc_tmp_path}/${Tumor}_${Normal}.${pos}.bam | grep "RG:Z:${Normal}" | awk -F'\t' '{print $5}' | ${datamash} median 1`

##############################################################
# 9、计算碱基质量
## datamash需要按列计算
## Tumor alt reads
base_quality=""
for base in `cat ${vcf_qc_tmp_path}/${Tumor}_${Normal}.mpileup.${pos}.list | awk -F'\t' '{print $6}' | sed 's/./&\n/g' | awk NF`
do
quality_tmp=`printf "%d\n" "'$base"`
quality=`expr ${quality_tmp} - 33`
base_quality=${base_quality}" "${quality}
done
median_quality_Tumor_alt=`echo ${base_quality} | tr ' ' '\n' | ${datamash} median 1`


## Tumor ref reads
${samtools} mpileup --min-BQ 0 --min-MQ 0 -B -d10000 --count-orphans -sf ${ref_fasta} --output-BP -r ${pos} \
${vcf_qc_tmp_path}/${Tumor}_${Normal}.${pos}.refreads.tumor.bam \
> ${vcf_qc_tmp_path}/${Tumor}_${Normal}.mpileup.${pos}.refreads.list 

base_quality=""
for base in `cat ${vcf_qc_tmp_path}/${Tumor}_${Normal}.mpileup.${pos}.refreads.list  | awk -F'\t' '{print $6}' | sed 's/./&\n/g' | awk NF`
do
quality_tmp=`printf "%d\n" "'$base"`
quality=`expr ${quality_tmp} - 33`
base_quality=${base_quality}" "${quality}
done
median_quality_Tumor_ref=`echo ${base_quality} | tr ' ' '\n' | ${datamash} median 1`


## Normal ref reads
${samtools} mpileup --min-BQ 0 --min-MQ 0 -B -d10000 --count-orphans -sf ${ref_fasta} --output-BP -r ${pos} \
${vcf_qc_tmp_path}/${Tumor}_${Normal}.${pos}.refreads.normal.bam \
> ${vcf_qc_tmp_path}/${Tumor}_${Normal}.mpileup.${pos}.refreads.normal.list 

base_quality=""
for base in `cat ${vcf_qc_tmp_path}/${Tumor}_${Normal}.mpileup.${pos}.refreads.normal.list  | awk -F'\t' '{print $6}' | sed 's/./&\n/g' | awk NF`
do
quality_tmp=`printf "%d\n" "'$base"`
quality=`expr ${quality_tmp} - 33`
base_quality=${base_quality}" "${quality}
done
median_quality_Normal_ref=`echo ${base_quality} | tr ' ' '\n' | ${datamash} median 1`



##############################################################
# 10、zeroproportion
## 比对质量为0的reads的数量
zeroproportion_Tumor=`${samtools} view ${vcf_qc_tmp_path}/${Tumor}_${Normal}.${pos}.bam | grep "RG:Z:${Tumor}" | awk -F'\t' '{if($5==0)print}' | wc -l`
zeroproportion_Normal=`${samtools} view ${vcf_qc_tmp_path}/${Tumor}_${Normal}.${pos}.bam | grep "RG:Z:${Normal}" | awk -F'\t' '{if($5==0)print}' | wc -l`

## 计算比例
zeroproportion_Tumor=`echo "" | awk '{print zeroproportion_Tumor/tumor_counts}' zeroproportion_Tumor=${zeroproportion_Tumor} tumor_counts=${tumor_counts} `
zeroproportion_Normal=`echo "" | awk '{print zeroproportion_Normal/tumor_counts}' zeroproportion_Normal=${zeroproportion_Normal} tumor_counts=${tumor_counts} `

##############################################################
# 11、strandbiasprop
# a,c represent the forward and reverse strand allele counts of REF reads and 
# b,d represent the forward and reverse strand allele counts of ALT reads. 
# |b/(a+b)-d/(c+d)|/((b+d)/(a+b+c+d))
strandbiasprop=`echo "" | awk '
{ 
strandbiasprop= ((b/depth_read1) - (d/depth_read2)) / (b+d)/(depth_read1 + depth_read2) ;
if(strandbiasprop > 0){strandbiasprop=strandbiasprop}
if(strandbiasprop < 0){strandbiasprop=-strandbiasprop}
print strandbiasprop
}' \
b=${Tumor_read1_alt_counts} d=${Tumor_read2_alt_counts} depth_read1=${Tumor_read1_counts} depth_read2=${Tumor_read2_counts}`

## 如果所有reads均在一个链上
if [ -z ${strandbiasprop} ]
then
strandbiasprop=1
fi

# strandbiassimple 
# min(forward/(forward+reverse),reverse/(reverse+forward)) 
strandbiassimple=`echo "" | awk '
{ 
value1=( depth_read1/(depth_read1 + depth_read2));
value2=( depth_read2/(depth_read1 + depth_read2));
if(value1 > value2){strandbiassimple=value2};
if(value1 <= value2){strandbiassimple=value1};
print strandbiassimple
}' \
depth_read1=${Tumor_read1_counts} depth_read2=${Tumor_read2_counts}`

##############################################################
## 12、突变counts数目在tumor和normal的比例
alt_counts_normal=`${samtools} view ${vcf_qc_tmp_path}/${Tumor}_${Normal}.${pos}.altreads.normal.bam | wc -l`
alt_counts_tumor=`${samtools} view ${vcf_qc_tmp_path}/${Tumor}_${Normal}.${pos}.altreads.tumor.bam | wc -l`

alt_vaf_normal=`echo "" | awk '{print alt_counts_normal/normal_counts}' alt_counts_normal=${alt_counts_normal} normal_counts=${normal_counts} `

##############################################################
## 总结输出
echo "${line},${tumor_counts},${normal_counts},${alt_counts_tumor},${alt_counts_normal},${alt_vaf_normal},${DistanceToAlignmentEndMedian},${DistanceToAlignmentEndMAD},${minmapqualitydifference},${medianmapquality_Tumor},${medianmapquality_altTumor},${medianmapquality_Normal},${median_quality_Tumor_alt},${median_quality_Tumor_ref},${median_quality_Normal_ref},${zeroproportion_Tumor},${zeroproportion_Normal},${strandbiasprop},${strandbiassimple}" \
>> ${vcf_qc_path}/${Tumor}_${Normal}_MutQc.csv
##############################################################
echo "Finish ${line} "
fi
fi




#### https://www.jianshu.com/p/ff6187c97155
## CIGAR标识符表示read中每个碱基的比对情况，主要有以下标识符
# M: alignment match (can be a sequence match or mismatch)
# read上的碱基与参考序列“RNAME”完全匹配，碱基一一对应，包括了正确匹配与错误匹配
# I: insertion to the reference
# read上的碱基相对于参考序列“RNAME”有插入现象（如下）：
# D: deletion from the reference
# read上的碱基相对于参考序列“RNAME”有删除现象（如下）：
# N: skipped region from the reference
# read上的碱基相对于参考序列“RNAME”存在连续没有比对上的空缺，这些空缺用N来表示，跟“D”相似但远比“D”缺失的更多，这种比对类型也叫“Spliced alignment”
# S: soft clipping (clipped sequences present in SEQ)
# H: hard clipping (clipped sequences NOT present in SEQ)
# read的开头或者结尾部分没有比对到参考序列"RNAME”上，但这部分未比对上的连续序列仍保留在sam文件的该read序列中，用“S”来表示；如果未保留，则用“H”表示，也即“hard cliping”（如下所示，也可同图2中r003的比对CIGAR中看出）
# P: padding (silent deletion from padded reference)
# 多条read比对到参考序列的同一位置时，如果不同read单独同该参考序列比对时，参考序列的情况也不同，
# =：sequence match 正确匹配
# X：sequence mismatch 错误匹配



## https://www.rdm.ox.ac.uk/files/research/lunter-group/stampyreadme.txt
# Multiple mapping locations
#   Normally, when Stampy identifies several equally good mapping 
# locations for a read or read pair, it reports one of these at random
#  (and assigns the choice a low mapping quality, of 3 or less).

#  Alternatively Stampy can report a limited number of alternative
#  mapping locations if you set the --xa-max option to a nonzero
#  value.  The alternative locations are reported in a single XA:Z:...
#  tag, using the format that is used by BWA.

#  Normally, no alternative mappings are reported for discordant read
#  pairs. However, if you set --xa-max-discordant to a nonzero value, 
#  these will be reported as well.


## gatk-tools的java脚本
## https://www.demo2s.com/java/apache-commons-descriptivestatistics-getelement-int-index.html